Average ground-state energy of finite Fermi systems 
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Semiclassical theories like the Thomas-Fermi and Wigner-Kirkwood methods give a good descrip- 
tion of the smooth average part of the total energy of a Fermi gas in some external potential when 
the chemical potential is varied. However, in systems with a fixed number of particles N, these 
methods overbind the actual average of the quantum energy as N is varied. We describe a theory 
that accounts for this effect. Numerical illustrations are discussed for fermions trapped in a har- 
monic oscillator potential and in a hard wall cavity, and for self-consistent calculations of atomic 
nuclei. In the latter case, the influence of deformations on the average behavior of the energy is also 
considered. 
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I. INTRODUCTION 

A basic problem in the physics of finite fermion sys- 
tems such as, e.g., atoms, nuclei, helium clusters, metal 
clusters, or semiconductor quantum dots, is the determi- 
nation of the ground-state energy E. A standard decom- 
position, deeply rooted in the connection of classical and 
quantum physics, is to write E as the sum of an average 
energy E and a fluctuating part E 0, 0, Q : 



E(N) = E{N) + E(N) 



(1) 



The largest contribution, E, is a smooth function of the 
number N of fermions. The shell correction E has a 
pure quantal origin and displays, instead, an oscillatory 
behavior as a function of N. 

Equation underlies the usefulness of the so-called 
mass formulae, like the liquid drop model for nuclei or for 
metal clusters, of which the oldest example is the well- 
known Bethe-Von Weizsacker mass formula for the bind- 
ing energy of nuclei. The decomposition Q is also at the 
basis of semiclassical and statistical techniques that are 
used to investigate how the properties of global charac- 
ter of fermion systems vary with the particle number N. 
Such is the case for instance of the celebrated Thomas- 
Fermi and Wigner-Kirkwood theories [HQ. These meth- 
ods often provide deep physical insights which may be 
otherwise obscured behind a full quantum calculation. 

It is recognized, however, that the semiclassical calcu- 
lations of E(N) for fermion systems in either external 
potentials or self-consistent mean fields show systematic 
deviations with respect to the actual average of the ex- 
act quantum results 0, 0, 0, IE IE 13 ■ F° r example, 
in spherically symmetric calculations one finds that, as 
a function of the number TV of particles, the difference 
E(N)-E(N) between the (fluctuating) exact value E(N) 
and the (smooth) semiclassical average E(N) does not 
oscillate around zero. In general, for fermions in a fixed 



external potential, semiclassical calculations overbind the 
true average of the quantum energy. One of our purposes 
in the present work is to explain the origin of this effect, 
and to derive an explicit formula that allows to compute 
the correct average behavior of E(N) in fermion systems. 
Related studies are the works of Refs. jSjEl where a par- 
ticle number conserving shell correction method has been 
pursued. 

Additional contributions to the average part of the 
ground-state energy come in fact from a careful analysis 
of the oscillatory term E(N). Because this fluctuating 
function is evaluated at discrete values of the chemical 
potential (that correspond to integer values of the par- 
ticle number), its average value is generically non-zero 
and therefore contributes to the average part of E(N). 
This phenomenon is related to the different physical de- 
scriptions of quantum mechanical systems obtained from 
different thermodynamic ensembles, the grand canonical 
and the canonical in the present context. This subtle 
topic has played, in recent years, a crucial role in under- 
standing the physics of, e.g., persistent currents in meso- 
scopic metallic rings 111. or in trapped Bose-Einstein 
condensates [T^ |. 

Our results are illustrated with two schematic models. 
Namely, we study the average of E(N) for fermions in a 
harmonic oscillator (HO) potential, via the semiclassical 
Wigner-Kirkwood (WK) theory jHj, and for fermions in 
a spherical cavity with sharp boundaries, via the Weyl ex- 
pansion |l4j| . In the former case, analytical expressions 
are available. Finally, and in contrast to the previous 
examples where the confining potential is fixed, we con- 
sider the influence of deformations and self-consistency 
on the average behavior of E(N), as well as other re- 
lated topics. We find that for self-consistent potentials 
with deformation degrees of freedom the behavior of the 
average energy is qualitatively different. 
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II. SMOOTH BEHAVIOR: GRAND 
CANONICAL VERSUS CANONICAL 
ENSEMBLES 

The usual computation of the different terms in Eq. 
is as follows. The single-particle level density g(e) = 
Tr[5(e — H)] of a quantum fermion system can be ex- 
pressed as 0, 0, E3 
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9(e) 



dfe(r,P) 
de 



dpdf 



(2) 



in terms of the phase space Wigner function f e (f,p). We 
have included a factor 2 to account for spin degeneracy. 
Then, for a set of fermions in a potential well filled up 
to an energy /i, the number of states (accumulated level 
density) and the ground-state energy are obtained from 
g(e) through 



g{e) de , E(p) 



eg{e)de . (3) 



Inserting the Wigner-Kirkwood expansion of the Wigner 
function f e (r,p) in powers of h in Eq. @ produces a 
smooth function g(e), where the leading order gives rise 
to the Thomas-Fermi term. T his p rocedure is well doc- 
umented in the literature P, 0, [H, 03 LL3 ■ Inserting 
the latter series for g(e) « g(e) into Eq. yields the 
semiclassical h expansions for Af(fi) and E(n). Alterna- 
tively, for a Fermi gas contained in a hard wall cavity, one 
inserts in Eqs. 10 the corresponding Weyl expansion EH 
of the average single-particle level density g(e). In both 
cases, Eqs. (0 produce a series in decreasing powers of 
whose coefficients depend on the shape of the potential. 

These expressions provide in general accurate descrip- 
tions of the average behavior of g(e), 7V(/i), and E(fi). 
For instance, for an isotropic three dimensional HO po- 
tential of frequency ui one obtains the well-known WK 
expressions H S El E El 13 
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The last term in Eq. Q contains the derivative of the 
delta function 5(e). This term and the last term in Eq. 
© stem from the corrections of order ft 4 to g and E, 
respectively. In the HO potential the h A contribution to 
AT vanishes. 

Figure ^ displays the comparison between the exact 
quantum mechanical quantities and the smooth approxi- 
mations (JSJ)-©. The upper panel shows the accumulated 
level density Af(n) for a set of fermions in a spherical HO 
potential, as a function of n/hco. The quantum result ex- 
hibits discontinuities at each major shell (N — 2, 8, 20, 
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FIG. 1: Accumulated level density (upper panel) and total 
energy (lower panel) with spin-degeneracy 2 for a spherical 
HO potential as a function of chemical potential /j,. Staircase, 
solid, and dashed lines correspond to the quantum, semiclas- 
sical WK, and shell correction (quantum minus semiclassical) 
values, respectively. 



40, 70, 112 in the present case) and is represented by a 
staircase function which fluctuates around the smooth 
WK curve provided by Eq. (J5J. The oscillatory part 
of Af(n) (dashed curve) contains the fluctuations due to 
shell effects. They are seen to oscillate around zero, with 
a vanishing net average, as /i is varied. The lower panel of 
Fig. displays the ground state energy E(p)/tuv for the 
same potential E3 . Again, the smooth WK curve excel- 
lently averages the quantum result and the shell energy 
fluctuates around zero. 

The fact that the average behavior of the remaining 
shell corrections is zero for E(n) and Af(fi) can be ex- 
plicitly checked. The general semiclassical theory ex- 
presses the fluctuating parts E(fj,) = E(ji) — E(fj.) and 
Af((i) = Af((i) — M(fj) as sums over the classical periodic 
orbits of the system at energy /j, 0, H3, 0, ^ . Each term 
in E(n) and Af(fi) is an oscillatory function of the chem- 
ical potential (through the action of the corresponding- 
orbit), whose average over a chemical potential window 
is zero. In the particular case of the HO potential the 
semiclassical approximation turns out to be exact (see 
e.g. 0,H3), and takes the form 
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where M r = 2irM, fi r = M r /i/Hu}, and Af(/J.) and E(/j,) 
are given by Eqs. JSJ and ©, respectively. These ex- 
pressions illustrate explicitly that the average of the fluc- 
tuating parts E([i) and Af(p) over a chemical potential 
window is zero, {E{^))^ = (/^(fi))^ = 0. 

In real physical Fermi gases with a well defined number 
of particles the various quantities, like masses or many- 
body level densities, are not studied as a function of the 
chemical potential /j, but rather as a function of the par- 
ticle number TV. For instance, the ground-state energy 
E(N) of the system consists in the sum of the single- 
particle energies of the N lowest single-particle states 
(taking into account spin-degeneracy). Thus Eqs. 
and (JSJ are related to the grand canonical ensemble. The 
qualitative behavior of the function E(N) as a function 
of N is in general quite different from the behavior of 
EM. 

Based on e.g. the Wigner-Kirkwood method, the usual 
way to calculate the function E(N) is as follows. Hav- 
ing determined, in that approximation, the energy E(p) 
and the accumulated level density A/"(/i), one first fixes 
the chemical potential (or Fermi energy) in terms of the 
particle-number N by inverting the function 

Af(p) = N^fl N = fl(N) . (9) 

To be consistent in the notation, we use fl N to denote 
the value of the chemical potential for a given N deter- 
mined from the WK (or Weyl) approximation, to stress 
that, in this approximation, is computed by inverting 
the smooth part AT and not the exact counting function 
Af. Finally, one obtains the smooth Wigner-Kirkwood or 
Weyl term E(N) replacing fi by p, N in E(n), 

E(N) ~ E(p, N ) = / eg{e) de . (10) 
Jo 

For example, applying this procedure to the isotropic HO 
potential, from the leading terms of Eqs. JSJ and © one 
straightforwardly obtains 

E(N) = \{2,N)^nio , (11) 

which is the leading-order Thomas- Fermi result. It shows 
that in a HO the leading dependence of the average en- 
ergy per particle, in units of huj, is oc TV 1 / 3 . 

The full Wigner-Kirkwood function E(N)/N com- 
puted for the HO potential including up to the fourth- 
order contributions in H is plotted in Fig. [3 (dashed line) 
as a function of the particle number N, in units of huj. It 
is compared to the exact quantum result (solid line). To 
better visualize the quantum oscillations with changing 
N, we have subtracted the dominant TV 1 / 3 dependence 
[recall Eq. 1)11(1 ] from both the quantum and the WK 
curves. The upper panel of Fig. |21 displays the results 
for the isotropic HO potential. The lower panel is for a 
strongly deformed potential and it will be discussed later 
on. Focusing on the isotropic HO, one sees that, as ex- 
pected, the general trend of the smooth WK result turns 




N 

FIG. 2: Upper panel: quantum and WK values of the en- 
ergy per particle for a spherical HO potential as a function 
of the number of particles, in units of hw. The leading N- 
dependence given by Eq. 1111 is subtracted from both curves. 
Lower panel: the same as in the upper panel but for a strongly 
triaxially deformed HO potential. Notice that the WK curves 
are different in the spherical and deformed cases. 



out to be quite correct in comparison with the global par- 
ticle number dependence of the quantum energies. There 
is, however, a systematic deviation in the sense that the 
WK curve does not pass as a function of N through the 
average of the quantum values. This is clearly seen from 
the large asymmetry of the shaded regions above and be- 
low the WK curve in the upper panel of Fig. [5] One 
notices that the WK result overbinds with respect to the 
true average of the quantum values when N is varied in 
the spherical symmetry. The same situation prevails in 
other problems of atomic and nuclear physics as well as in 
self-consistent mean field calculations 0, IE EL IE IE 13 • 



III. CONTRIBUTION OF THE OSCILLATORY 
CORRECTIONS 

The previous results show that the function E(N) does 
not describe appropriately the average behavior of E(N). 
We now discuss the origin of the discrepancy, and the way 
to correct it. 

In systems with a well defined number of particles the 
chemical potential fi takes discrete values. These values 
do not occur at random. For instance, for an even num- 
ber of particles and non-degenerate single-particle states, 
a standard rule is to locate the chemical potential half- 
way between the last occupied and the first unoccupied 
single-particle states. Fixing a particular rule to deter- 
mine the chemical potential at a given number of particles 
introduces a bias in the sampling of the values of \x (with 
respect to a uniform, random distribution of /i). Due to 
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this bias, when the oscillatory part of the energy E(ijl) 
is evaluated over the set of discrete points it produces, 
generically, a function whose average is different from 
zero. To compute that average we proceed as follows. 

The decomposition of the single-particle level density 
into a smooth part and a fluctuating part, 

3(e) =fl(e)+fl(e) , 

where g(e) is the WK (or Weyl) smooth part and g(e) is 
given by the sum over periodic orbits mentioned above 
1], induces a corresponding decomposition of the inte- 
grated density [cf Eq. J3J)] : 

7V(m) =#(/*) +#0*) . (12) 

For a given number of particles N, the chemical po- 
tential is defined by inversion of the exact accumulated 
level density 

7V(m) = N fi N = fx(N) . (13) 

As the particle number N increases, it is natural to de- 
compose the chemical potential into smooth and fluctu- 
ating parts: 

The average part p, N satisfies Eq. ©. Assuming that 
/ijv <C p, N , a Taylor expansion of the smooth part in 
powers of JL N around /i = p, N in Eq. (|12fl (no Taylor 
expansion is allowed for the fluctuating term Af, because 
it is not a regular function) yields, to lowest order, 

V>n = -\ A^>«) , (14) 
9 

where we denote 

g = g(p, N ) . (15) 

Similarly, the energy may be decomposed as 

Efa) = E(n)+E{n) . (16) 

In a system with a well-defined number of particles, the 
smooth part E(N) of the exact function E(N) = E(fj, N ) 
was defined in CDJ, E(N) = E(p, N ) = f£ N e g(e) de. The 
fluctuating part is thus defined as 

E(N) = E(N) - E(N) . (17) 

In order to compute E(N), and in particular its aver- 
age over some particle-number window AN around N, it 
is convenient to express the energy in terms of the grand 
potential f2 = — f ll J\f(e)de using the thermodynamic re- 
lation 

E(ji) = SlQi) + ntfQi) . 



Recalling the definition of fi N and fi N [Eqs. © and (|13|) ]. 
Eq. H17|) may be written 

E(N)=Sl(fi N )-n(p N ) + Jl N N , 

where fi(/2 N ) = — f 1 *™ M(e)de. Decomposing Q(fi N ) 
into its average and fluctuating parts, expanding Q(fj, N ) 
around p, N to second order in J1 N , and using the thermo- 
dynamic relations 

— — =-M(ti N ) and = -g , 

we get 

E(N)=ti( f i N )-gji 2 N /2 . 
Using Eq. Q14|). this takes the form 

E(N) = h(fx N ) - -L Kf 2 {^ N ) + OQll) . (18) 
2 9 

Equation l|18|) connects the fluctuations of the grand po- 
tential (grand-canonical ensemble) to those of the energy 
at a fixed number of particles (canonical ensemble). This 
connection, to lowest order, has been exploited in recent 
years to analyze nuclear-mass fluctuations [241 l2fij . 

One may be tempted to think that the average of E(N) 
over some particle-number window AiV around N, de- 
noted (E(N)) N , is proportional, from Eq. ljT8)l. to the 
variance (Af 2 (fJ. N )) N of N and that the average of 
E(N), 

(E{N)) N = E(N) + (E(N)) N , (19) 

is thus lowered with respect to E(N) (due to the minus 
sign in front of Af 2 (/J, N ) in Eq.JTS}). However this is 
wrong because, for the same reasons as for E(N), the 
average (fl(ii N )) N ^ 0. A detailed calculation (cf the 
Appendix) shows that 

= r <AA 2 (^))„ + I 9 (4>« " , (20) 
9 8 6 g 

where (s 2 )« is the variance of the spacing s N = e N+1 — e N 
between two consecutive single-particle levels around /j, N . 
Taking the average with respect to the discrete points fi N 
in Ea. l|18[l . using Eq. (|20|) for the average of (l(fi N ) and 
expressing, for convenience, the average (Af 2 (fi N )) N over 
the discrete points /i N in terms of the average (A r2 (/i)) A1 
over the continuous variable fi around fj, N (cf the Ap- 
pendix), 

{Af 2 M) N = (Af 2 ^)), + \-\f (4>« , (2i) 

we get 

(E(N)) N = ± (Af 2 (fi))n + 0(Jll) . (22) 
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FIG. 3: Quantum and WK values of the energy per particle 
in a triaxially deformed HO potential. Spin degeneracy is 
included. 



The final expression for the average value of the energy in 
a system conserving the number of particles is, according 
to Eqs. and J22J), 



1 



(E(N)) N — E(N) + — = (Af (/i)) 
2 g 



i 



12 g 



+ . (23) 



It follows from Eq. (|2"3l) that, with respect to the WK 
or Weyl smooth terms, the true average of the energy as 
a function of N is increased by a term proportional to 
the variance of the accumulated level density. Equation 
(I23[) contains all the relevant information on the average, 
and allows to understand the numerical results presented 
above. Before making a quantitative comparison, we first 
discuss the general aspects involved in that equation. 

Equation (|20|l is demonstrated in the Appendix for 
a system without degeneracies (intrinsic and/or due to 
spin). However, it is easy to see that it is also valid in 
the presence of degeneracies. This is because the ther- 
modynamic quantities we are considering are continuous 
variables of a given set of external parameters A. As- 
sume that for some A = An. degeneracies occur, and that 
for slightly different values A ^ A all the degeneracies 
are lifted (for instance, some of the components of A may 
be associated to a shape deformation, with A = An the 
spherical case, and another component may be a mag- 
netic field that lifts the spin degeneracy, with A = An = 
no magnetic field). Then for A ^ An Eq. Q20(l is valid. 
One can therefore consider the case with degeneracies as 
the limit A — + An and, by continuity, Eq. I|2U|) remains 
valid. 

The variance (A/" 2 (^)) M in Eq. depends on the sys- 
tem under consideration. However, its general properties 
can easily be determined. In systems where the typi- 
cal size of the fluctuations is important, then the shift 
of the true average (E(N)) N with respect to E(N) will 
also be important. On the other hand, in systems with 
small fluctuations, (N 2 (n))^/2g will be small, and the 
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FIG. 4: Upper panel: normalized fluctuating part E(N) = 
E(N) — E(N) as a function of N for a spherical cavity. 
Eo = h 2 /2mj-Q where rn is the radius of the sphere and m 
is the mass of fermions. Lower panel: average {E(N)) N of 
the fluctuating function in the upper panel (full line) com- 
pared to the theoretical prediction, Eq. 122H (dashed line). 



term E(N) will give not only a good approximation to 
(E(N)) N , but also to E(N) as well (since fluctuations are 
small). In general, the more regular and/or symmetric a 
system is, the greater the fluctuations are, and the larger 
the correction (A/" 2 (/x)) )Jl /2g will be. As the regularities 
or symmetries are broken, the typical size of the fluctu- 
ations diminishes, and E(N) will be well approximated 
by E(N). This point is illustrated in Fig. |2 where the 
upper panel shows E{N)/N for the isotropic HO, where 
large fluctuations are observed and clear deviations of 
the average with respect to E(N) are found. In contrast, 
the lower panel shows a strongly deformed HO, with fre- 
quencies uj x /u> — 0.460, LOy/oj — 1.111, and lu z /lu = 1.954 
{uj x L0 y u) z — u! 3 ), where small fluctuations are observed, 
and as a consequence good agreement between E(N) /N 
and E{N)/N is found. Another manifestation of the 
same phenomenon is provided in Fig. where E(N) /N 
is shown for N = 92 fermions (with spin degeneracy) 
in a triaxially deformed HO potential as a function of 
the deformation parameter d, where uj x /uj = <5 _1 / 2 /cr 1 ' /3 , 
ujy/uj = 5 1/2 /o- 1/3 , and lj z /uj = a 2/3 , with a = 1 + d\/3 
and <5 = 1 + \d\y/2. We sec that for most deformations 
(mid-shell configurations) the quantum and the smooth 
WK values practically agree, up to small fluctuations. 
Large deviations are observed, instead, when sphericity 
is approached, and for other special deformations, e.g., 
for d ~ 0.65 where the frequency ratio lo x :lo v :oj z is close 
to 1:2:3 (when the three frequencies ui x , uj y , and tu z are 
integer ratios, the energy levels of the HO are degenerate 
and the classical trajectories of the Hamiltonian become 
closed periodic orbits |l(). 

We have made a quantitative check of Eqs. 1221) and 
(|23|l for the case of a Fermi gas in a spherical cavity. 
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FIG. 5: Upper panel: fluctuating part E(N) = E(N)-E(N) 
as a function of N for a 3D isotropic HO. Lower panel: average 
{E(N)) N of the fluctuating function in the upper panel (full 
line) compared to the theoretical prediction, Eq. 1261 (dashed 
line) . 

The upper panel of Fig. 0] represents the fluctuating 
part E(N) as a function of N, defined in Eq. ljT7)l. A 
clear structure organized in shells (rapid oscillations) and 
supershells (long-range modulation of the rapid oscilla- 
tions) is observed. The lower panel shows a compari- 
son between the average (E(N)) N calculated numerically 
from the upper panel of that figure and the result ob- 
tained from Eq. (1221 as a function of N. The average 
shows a nontrivial dependence with the particle number 
(which reflects, to a large extent, the supershell struc- 
ture), that is very well reproduced by theory. 

In the case of the spherical HO, it is possible to get eas- 
ily an analytical expression for (E(N)) N . The function 
Af(/i) is given by the second term in the right hand side 
of Eq. |0| . Squaring it, the main contribution to A/" 2 (/i) 
comes from terms where both indices of the double sum 
are equal. Hence to leading order in we get: 

^)^(£) 1 f:(« ZM ) 2 - m 

v 7 M=l v 7 

(A/' 2 (/i)) Al is calculated by averaging over the rapidly fluc- 
tuating factors, given by the sine terms. This yields 

M=l v 7 

Since (A/" 2 (/i)) M is a smooth function, we replace \x by ft. 
Thus we can use the WK expression Eq. © to compute 
the dependence of the variance with the number of par- 
ticles. Using moreover Eq. 10} we finally get, to leading 
order in A, 

(E(N)) N = ^ith = 1(37V) 2 / 3 ^ . (26) 



A comparison with the numerical average of E(N), ob- 
tained from an isotropic 3D HO, is presented in Fig. [5] 
The result shows an excellent agreement; compared to 
the spherical cavity, a much simpler N dependence is ob- 
served, due to the absence of supershells. 

Based on general properties of the single-particle spec- 
trum, it is possible to estimate the typical size of the vari- 
ance (A/" 2 (/x)) M and of its A dependence for a large class 
of systems, and therefore to estimate (E(N)) N . The rele- 
vant classification relies on the type of classical dynamics 
associated to the confining potential. The two extreme 
cases that can be treated explicitly are fully regular and 
fully chaotic dynamics (the case of mixed dynamics is 
more subtle). Based on this classification, explicit re- 
sults for the typical size of (AT 2 (/x)) M were obtained in 
Ref. [13. 

IV. DEFORMED AND SELF-BOUND SYSTEMS 

Up to now we have considered fermion systems con- 
fined by an external potential. This may be applica- 
ble to e.g. quantum dot systems or magnetically trapped 
atomic gases where the self-consistent mean field part 
plays a minor role with respect to the external confin- 
ing potential. However, many relevant systems are self 
bound and then the mean field potential is essentially 
given by the solution of the self-consistent Hartree-Fock 
equations which are obtained by minimising the energy 
of the system. The mean field in these situations may 
turn out to be spherical, but in many cases rotational 
invariance will be broken and the mean field becomes de- 
formed. We will see that in these cases the results can 
show interesting differences with regard to the scenario 
found in the upper panel of Fig. [21 or in Fig. |3] for the 
harmonic oscillator, and in Fig. 0] for the hard wall cav- 
ity, where the potential was kept spherical. We want to 
investigate such cases now. 

First, to illustrate the situation, we again consider 
the HO potential. In contrast to the previous section, 
for each particle number we now minimise the ground- 
state energy of the quantum solution with respect to 
deformation, i.e., with respect to free variation of ui x , 
u) y , and uj z , under the constraint of volume conservation 
(u) x oj y u> z = ui 3 ). This must be done in carefully check- 
ing simultaneously the optimal choice of the occupan- 
cies n x ,n y ,n z . The semiclassical energies E(N) always 
have their absolute minimum at sphericity as given by 
Eq. ©• As a particular example of self-bound system, 
we consider the case of the atomic nuclei. We mimic 
the saturation properties of nuclear forces by including 
the standard particle-number dependence of the HO fre- 
quency hw = 41 A- 1 / 3 MeV @ with A = 2N (i.e., A here 
represents the mass number of a hypothetical uncharged 
nucleus with A protons and A neutrons). 

In Fig. we show the difference 5E between the 
fully minimised quantum energies and the corresponding 
isotropic semiclassical expression E{N) obtained from 
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FIG. 6: Difference between the minimised ground state en- 
ergy with respect to deformations and the spherical WK en- 
ergy E(N) for fermions in a HO potential (squares joined by 
full line). For comparison, the dotted curve shows the energy 
difference for an isotropic HO (same curve as in the upper 
part of Fig. |5J. The average of both curves is shown by a 
long-dashed line. The scaling hu> = 41(27V)" 1/3 MeV has 
been used in the calculations. 



Eqs. © and iJBJ. For comparison, we include in the 
same figure the fluctuating part E(N) for the spherical 
HO (same curve as in the upper part of Fig. [5}. Wc 
observe that with respect to the purely spherical case, 
the situation changes very much. Now, contrary to the 
spherical case, practically all values of SE from the min- 
imised quantum solutions are negative, meaning that the 
minimised quantum energies stay below the semiclassi- 
cal curve E(N). The minimised quantum energies coin- 
cide with the spherical ones only in a small neighborhood 
around the shell closures, whereas away from the latter 
the system is axially deformed or even, around the middle 
of the shells, a slight triaxiality can appear (in the case of 
axial symmetry, typical deformations show an axis ratio 
of 2:3). 

It seems natural that the deformed quantum ener- 
gies are more bound than the approximate energies ob- 
tained from the semiclassical theory, in spite of the fact 
that to our knowledge no upper bound theorem like the 
Rayleigh-Ritz principle exists for the semiclassical ap- 
proach. We wish to point out that in Fig. for most val- 
ues of N the system is actually rather well deformed and 
that, with the exception of a couple of particle numbers 
around closed shells, the energy differences are in most 
cases very close to the zero line, i.e., to the semiclassical 
WK values. This is consistent with the results obtained 
in the previous section, where it was shown that for de- 
formed systems where degeneracies are lifted the energy 
E(N) is expected to be well approximated by the WK 
theory. 

The magnitude of the difference SE of the minimised 
quantum solutions to the semiclassical values slightly in- 
creases with increasing particle number, as the average 



curve shows in Fig. H3 However, the magnitude of the 
same quantity SE divided by the particle number de- 
creases as a function of increasing N, and the minimised 
deformed quantum energies per particle are extremely 
close to the semiclassical ones. Notice the opposite trend 
in Fig. of the two average curves, with and without 
energy minimisation. In contrast to the latter case, for 
which an explicit formula for the average behavior was 
developed and successfully checked in the previous sec- 
tion, we do not yet have an equivalent result for a self- 
bound system. 

We are interested in checking whether this simplified 
harmonic oscillator scenario remains valid in realistic 
Hartree-Fock-type mean field calculations. In Ref. @ 
self-consistent calculations of the ground-state binding 
energy of atomic nuclei were carried out using the varia- 
tional Wigner-Kirkwood method [271 • The nuclear inter- 
action was described by the relativistic mean field (RMF) 
meson exchange model 28]. Quantum calculations for 
the RMF model are available in the literature. In partic- 
ular, a mass table of deformed (axially symmetric) quan- 
tum calculations for nuclei with an accurately calibrated 
RMF nuclear interaction is published in Ref. (2^. From 
this table we took for each value of the mass number A 
the most bound (in general deformed) isotope and traced 
E/A as a function of A [3(j ■ The quantum values to- 
gether with the RMF semiclassical results, computed fol- 
lowing Ref. fH, are shown in Fig.JTJfor nuclei with an even 
number of protons and neutrons. Most of the WK ener- 
gies lie on top of the deformed quantum energies on the 
scale of the figure. We plot in addition the RMF quan- 
tum values constrained to sphericity. The typical arch 
structure found in Fig. [21 for the spherical HO poten- 
tial is then recovered. These arches in nuclei take place 
between the so-called magic numbers, i.e., the proton or 
neutron numbers where effects analogous to the shell clo- 
sures of the HO or of the electron shells in atoms occur. 
The fact that for nuclei above iron E/A raises whereas 
in Fig. |2 it keeps decreasing is a trivial effect due to the 
Coulomb repulsion among protons in the atomic nucleus. 

In Fig. [S] we display for the self-consistent RMF the dif- 
ference SE between the quantum energies (that are, as 
mentioned, minimised with respect to deformation) and 
the corresponding semiclassical values (that attain their 
absolute minima at sphericity). For reference, also the 
values of the energy difference SE obtained from experi- 
mental data |3lj instead of the quantum results are dis- 
played. The similarity of Fig. |S| with Fig.[S]for the HO is 
striking. Systems with the largest deformations are again 
located mostly around mid shells. When the system ap- 
proaches the spherical shape, SE becomes increasingly 
negative and displays the downward peaks seen in Fig. El 
on reaching neutron or proton magic numbers. 

It is clear that in the self-consistent case as in the 
schematic case of the HO with optimized shapes con- 
sidered above, the average of SE as a function of particle 
number is, at least for the heavier systems, negative. In 
self-bound systems we again find that in between shells 
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FIG. 7: Energy per particle of atomic nuclei with an even 
number of protons and of neutrons along the periodic table, 
as obtained from self-consistent relativistic mean field calcu- 
lations. The deformed calculations are from Ref. |29(. For 
each value of the mass number A we plot the most bound 
isotope according to the tabulation of Ref. |2{| (excepting for 
A = 40 [U). 



FIG. 8: Energy difference between the deformed quantum 
solutions and the WK values of Fig. 7 for the self-consistent 
relativistic nuclear mean field. The result obtained by taking 
the difference of the experimental data [3l| to the calculated 
WK energies is also shown for the purpose of illustration. The 
location of the magic neutron (N) and proton (Z) numbers 
is indicated. 



the quantum energies are closer to the semiclassical WK 
values. In between shells the system deforms in search of 
minimum energy and avoids the large positive shell cor- 
rections to the energy that occur if a spherical shape is 
kept. As the deformation increases, symmetries are bro- 
ken and the amplitude of the shell corrections diminishes. 
This is in agreement with the basic ideas underlying Eq. 
(I23[l . that imply an energy E(N) which is well approxi- 
mated by the WK theory. 



V. CONCLUDING REMARKS 

In this work we have revisited the old problem of the 
semiclassical approach to finite fermion systems, based 
either on the Wigner-Kirkwood expansion with ft correc- 
tions, or the Weyl expansion. We have addressed the 
nature of one of the most elusive features of the the- 
ory, namely the systematic overbinding compared to the 
quantum average for fermions in fixed potentials. 

In the first part, we have shown that this discrepancy 
is due to the fact that these methods do not incorpo- 
rate appropriately the conservation of the particle num- 
ber. There is, generically, a contribution to the aver- 
age ground-state energy that comes from the fluctuating 
part, or shell contribution. We derived an explicit for- 
mula that takes into account that contribution, and have 
tested it for different fixed confining potentials. In all 
cases, a positive correction with respect to the semiclas- 
sical result is predicted [cf Eq. |J23J], whose magnitude 
depends on the size of the shell effects. When the con- 
fining potential has symmetries, the shell corrections are 
large, and important deviations between the exact quan- 
tum and the WK energies are observed, in agreement 



with our predictions. In contrast, when symmetries are 
broken shell effects are smaller, and the exact energies 
(and not only their average part) are better described by 
the WK theory. 

The description of the behavior of self-bound systems 
is more difficult and subtle, because at each particle num- 
ber N the energy is minimised, and hence the shape of the 
potential is a function of N. In this shell correc- 

tion which is nearly always negative with respect to the 
spherical WK result is observed for the HO potential. In 
between shells, when the system deforms and symmetries 
are broken, the value of the shell correction is smaller, 
and the energy is well approximated by the WK theory, 
in agreement with the general considerations that follow 
from Eq. (|23[l . Interestingly, all these features have been 
qualitatively confirmed by a more realistic model based 
on a mean-field self-consistent calculation of the ground- 
state energy of atomic nuclei. However, the problem of 
deriving an explicit formula for the average behavior of 
the ground-state energy of self-bound systems that cor- 
rectly takes into account the iV-dependence with defor- 
mation degrees of freedom is still open. 

We are indebted to R. K. Bhaduri and O. Bohigas 
for valuable informations. Work partially supported 
by IN2P3-CICYT. X.V. and M.C. acknowledge Grants 
No. FIS2005-03142 (MEC, Spain, and FEDER) and No. 
2005SGR-00343 (Generalitat de Catalunya). P. L. and 
J. R. acknowledge support from grants ACI Nanoscience 
201, ANR NT05-2-42103, ANR-05-Nano-008-02 and the 
IFRAF Institute. 
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APPENDIX 



We follow here Appendix B of Ref. |32j in order to 
prove Eqs. and l|2T|). 

Let us consider a single-particle spectrum £,-, with 



J 



1,2,... and < e 



j+i- 



The accumulated level den- 



sity N{\i) is discontinuous at each energy level. At the 
discontinuity, we assign to the "intermediate" value 
Af(e N ) — A - 1/2. For A > 1 and e N < fi< e N+u writ- 
ing A/"(//) = A/"(/Lt) + , making a Taylor expansion of 
M(fJ.) around e N , and using that Af(e N ) — A — 1/2, we 
may write 

=Af(e N ) + -- {fi-e N )g . (Al) 

Evaluating this relation just before /i = s N+1 , and taking 
into account the value of the function at e N+1 , we have 



Af(e N+1 ) = Af (£«) + 1 - s N g 



(A2) 



where 



is the level spacing (we neglect the dependence of g with 
energy). Taking the discrete average over A on both 
sides of Eq. (|A2|> . defined as 



N+AN/2 

— y 

AN ^ 

j—N — AN/2 



where A A is the number of levels in the window around 
the Ath level, and using that (A/"(ejv + i)) N = (Af(e N )) N , 
we obtain the (trivial) relation (s N ) N = 1/g. 

On the other hand, defining the continuous average 
over a window A/i of a function that depends on the 
chemical potential as 



N+AN/2 ,. 



e j+i 



j=N—AN/2 3 



where AA = gA.fi is the number of levels in the window, 
we have from Eq. IjAlfl 



1 

A/i 



mi*))* = 

N+AN/2 

E 

j=N-AN/2 

AN (,rrt ^ 
g ( {Kf{e N )s N ) N 



£ j + l r 



N{£j) + 1/2 - g(fi - £j) 



dfi 



1 

2i 



2 (4, 







(A3) 
by defini- 

2— 2, . (M) 

Squaring and computing the discrete average in both 

sides of Eq. i|A2|) it is possible to deduce that (Af(e N )) N = 
after using Eq. (|A4|) . 



The last equality follows because (A/"(/i)) M 
tion. From Eq. (|A3|) we deduce 

1 



9 , 



N )°N I N 



Similarly, squaring Eq. (|A1|) and computing in both 
sides the continuous average it is possible to relate the 
continuous variance of Af(fi) with discrete averages 

<A^)) M = g(rt 2 (e N ) SN ) N ~g 2 {£f(e N ) S 2 N ) N +^( S 3 N ) N ~ . 

(A5) 

Computing the discrete average of the third power of 
Eq. (|A2|) , the previous expression for the continuous vari- 
ance is considerably simplified and gives 



(A6) 



We are now interested in the statistics at the discrete 
points fi N = (e N+1 + ejv)/2. From Eq. (|A1|) we have 

N{n H )=N(e N ) + ^-^s N , (A7) 

from which it is easy to deduce that (Af(fJ, N )) N = 0. From 
the discrete average of the square of Eq. (|A7|I , using the 
result <|A4j) as well as Eq. (|A6|) , it follows that 



6 4 



(A8) 



We have demonstrated Eq. l(2*T|) . 

Let us now consider the grand potential for A ^ 1 and 
e N < fi < e N+1 . Since Q(fx) — fl{e N ) = — N{e) de, 
using Eq. (|A1|) and integrating we get 



= n(e N ) - \ Af{e N ) + - J (ji - e N ) + |(/x - e N ) 2 . 

(A9) 

Noting that 0(/z) is a continuous function, Eq. i|A9(l at 

/i = £„+! gives 

ti{e N+1 ) = fl(e N ) - (ti(e N ) + ^ s N + |«* . (A10) 

In a similar way as it was done for the accumulated 
level density, by integration of Eq. (|A9|1 (knowing that 
(f2(/ x)) M = 0), and discrete average of the product of 
Eq. {A2j) and Eq. l|A10)) . we deduce 



(n(e N )) N = -(Af 2 ^)),, 
From Eq. I|A9(1 . for fi = fi N , we have 



(All) 



n(/x w ) = fi(e w ) - (AA(e N ) + 2 J Y + | s « ■ (A12) 

The discrete average of this equation, together with 
Eqs. I|A4|I and HA11|) . finally leads to 



1\ s„ , g 2 



(^«)) N = ^(AA 2 ( Ai )), I -f(4>« • 
5 ° 



(A13) 



This equation, together with Eq. i|A8|l . implies Eq. J^UJl. 
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